%%Solving for the transitional dynamics of a 3 country model



@#define  n=3
@#for i in 1:n
var rho_chi@{i}   B@{i} gyy@{i} Y@{i} YT@{i} YNT@{i}  C@{i} gt@{i} gtT@{i} gtNT@{i} Hr@{i} gz@{i}  gp@{i} P@{i} PT@{i} PNT@{i} Pi@{i} PiT@{i} PiNT@{i} w@{i} Vinnov@{i}  Z@{i} T@{i} TNT@{i} Y_L@{i} gy@{i} gc@{i} ;

@#endfor

var EX1 EX2 IM1 IM2 Yw r ;

@#for i in 1:n
@#for j in 1:n
var     chivar@{i}@{j} ga@{i}@{j}  A@{i}@{j} pi_share@{i}@{j} pi_shareNT@{i}@{j}   V@{i}@{j} Vinnovator@{i}@{j} rp@{i}@{j} ;
varexo distexo@{i}@{j} distexoNT@{i}@{j} chiexo@{i}@{j} ;
@#endfor
@#endfor


@#for i in 1:n
@#for j in 1:n
//@#if j != i
var   drp@{i}@{j} Ha@{i}@{j}  J@{i}@{j} Jinnovator@{i}@{j}  eps@{i}@{j};
//@#endif
@#endfor
@#endfor



// Model Parameters 


parameters  mbar betar sig betaa bet eta gam_barg;

@#for i in 1:n
parameters ipr@{i} ipr_rho@{i}  Bbar@{i} lam@{i} L@{i}  gam@{i} al@{i} LNT@{i} LT@{i} A@{i};
@#endfor

@#for i in 1:n
@#for j in 1:n
parameters   tau@{i}@{j} dist@{i}@{j} distNT@{i}@{j} epsbar@{i}@{j};
@#endfor
@#endfor

load Parameters.mat;
load Initial_values.mat;
//load Q.mat;
//load epsbar.mat;

mbar=1.25;
betar=0.5223;//0.4614;//0.3432;//0.2882;//0.4285;//0.4597;//0.3155;//692;//0.3931;//0.4622;//0.4450;//0.3883;//0.4661;//0.5904;//0.5918;//0.5497;//0.5301;//0.5424;//0.5562;//0.5519;//0.4825;//0.3;//0.36;//0.3316;//0.3673;//0.4999;
bet=0.985;
sig=5;
betaa=betar;//0.5;
eta=0.001;
gam_barg = 1;//0.3010;

//chi11 = chi(1,1);
///chi12 = chi(1,2);
//chi13 = chi(1,3);

//chi21 = chi(2,1);
//chi22 = chi(2,2);
//chi23 = chi(2,3);

//chi31 = chi(3,1);
//chi32 = chi(3,2);
//chi33 = chi(3,3);

ipr1 = ipr(1);
ipr2 = ipr(2);
ipr3 = ipr(3);

ipr_rho1 = ipr_rho(1);
ipr_rho2 = ipr_rho(2);
ipr_rho3 = ipr_rho(3);

lam1=lam(1);
lam2=lam(2);
lam3=lam(3);

Bbar1=0.;//Bbar(1);
Bbar2=0;//Bbar(2);
Bbar3=-0.;//Bbar(3);

A1=Q(1);
A2=Q(2);
A3=Q(3);

al1=alpha(1);
al2=alpha(2);
al3=alpha(3);


L1=L(1);
L2=L(2);
L3=L(3);


LT1=LT(1);
LT2=LT(2);
LT3=LT(3);

LNT1=LNT(1);
LNT2=LNT(2);
LNT3=LNT(3);

epsbar11=epsbar(1,1);
epsbar12=epsbar(1,2);
epsbar13=epsbar(1,3);

epsbar21=epsbar(2,1);
epsbar22=epsbar(2,2);
epsbar23=epsbar(2,3);

epsbar31=epsbar(3,1);
epsbar32=epsbar(3,2);
epsbar33=epsbar(3,3);


dist11=distT(1,1);
dist12=distT(1,2);
dist13=distT(1,3);

dist21=distT(2,1);
dist22=distT(2,2);
dist23=distT(2,3);

dist31=distT(3,1);
dist32=distT(3,2);
dist33=distT(3,3);



tau11=tau(1,1);
tau12=tau(1,2);
tau13=tau(1,3);

tau21=tau(2,1);
tau22=tau(2,2);
tau23=tau(2,3);

tau31=tau(3,1);
tau32=tau(3,2);
tau33=tau(3,3);




distNT11=distNT(1,1);
distNT12=distNT(1,2);
distNT13=distNT(1,3);

distNT21=distNT(2,1);
distNT22=distNT(2,2);
distNT23=distNT(2,3);

distNT31=distNT(3,1);
distNT32=distNT(3,2);
distNT33=distNT(3,3);





gam1=lamNT(1);
gam2=lamNT(2);
gam3=lamNT(3);


model;



//Resource Constraint

@#for i in 1:n
exp(Y@{i})=exp(C@{i})+exp(Hr@{i})
@#for j in 1:n
+exp(Ha@{i}@{j})
@#endfor
;
@#endfor


//Euler Equation

@#for i in 1:n
rho_chi@{i} = 0.25*(ipr_rho@{i})^gam_barg;
1+eta*(B@{i}*exp(Y@{i})-Bbar@{i})= bet*(1+gc@{i})^(-1)*exp(r);
@#endfor

B1*exp(Y1)+B2*exp(Y2)+B3*exp(Y3)=0;


//Prices

@#for i in 1:n
exp(PT@{i})^(1-sig)=
@#for j in 1:n
+((A@{j}))^(sig-1)*(exp(T@{j}(-1)))*(mbar*exp(w@{j})*dist@{i}@{j}*(1+tau@{i}@{j}*(1-distexo@{i}@{j})))^(1-sig)
@#endfor
;
@#endfor


@#for i in 1:n
exp(PNT@{i})^(1-sig)=
@#for j in 1:n
+((A@{j}))^(sig-1)*exp(TNT@{j}(-1))*(mbar*exp(w@{j})*distNT@{i}@{j}*(1-distexoNT@{i}@{j}))^(1-sig)
@#endfor
;
@#endfor



@#for i in 1:n
exp(P@{i})=(exp(PT@{i})/al@{i})^al@{i}*(exp(PNT@{i})/(1-al@{i}))^(1-al@{i});
@#endfor

//Trade shares


@#for i in 1:n
@#for j in 1:n
 exp(pi_share@{i}@{j})=((A@{j}))^(sig-1)*exp(T@{j}(-1))*(dist@{i}@{j}*(1+tau@{i}@{j}*(1-distexo@{i}@{j})))^(1-sig)*(exp(w@{j})*mbar)^(1-sig)/exp(PT@{i})^(1-sig);
 exp(pi_shareNT@{i}@{j})=((A@{j}))^(sig-1)*exp(TNT@{j}(-1))*(distNT@{i}@{j}*(1-distexoNT@{i}@{j}))^(1-sig)*(exp(w@{j})*mbar)^(1-sig)/exp(PNT@{i})^(1-sig);
@#endfor
@#endfor


//Labor market clearing condition

@#for i in 1:n
mbar*exp(w@{i})*LT@{i}=
@#for j in 1:n
+exp(pi_share@{j}@{i})*exp(YT@{j})*exp(P@{j})/(1+tau@{j}@{i}*(1-distexo@{j}@{i}))

@#endfor
;
@#endfor

@#for i in 1:n
mbar*exp(w@{i})*LNT@{i}=
@#for j in 1:n
+exp(pi_shareNT@{j}@{i})*exp(YNT@{j})*exp(P@{j})
@#endfor
;
@#endfor

// Output

@#for i in 1:n
exp(Y@{i})=exp(YT@{i});
@#endfor

exp(Yw)=exp(P@{1})*exp(Y@{1})+exp(P@{2})*exp(Y@{2})+exp(P@{3})*exp(Y@{3});

//Profits

@#for i in 1:n
exp(Pi@{i})= (1/(sig-1))*exp(w@{i})*L@{i}; 
exp(PiT@{i})= (1/(sig-1))*exp(w@{i})*LT@{i}; 
exp(PiNT@{i})= (1/(sig-1))*exp(w@{i})*LNT@{i}; 
@#endfor

//Value of adopted technologies

@#for i in 1:n
@#for j in 1:n
exp(V@{i}@{j})=(1-exp(chivar@{i}@{j}))*exp(PiT@{i})+(1/exp(r))*exp(V@{i}@{j}(+1))*((1+gt3)^(1/(sig-1))/(1+(gt@{i})))*(exp(P@{i})/exp(P@{i}(+1)));
exp(Vinnovator@{i}@{j})=exp(chivar@{i}@{j})*exp(PiT@{i})+(1/exp(r))*exp(Vinnovator@{i}@{j}(+1))*((1+gt3)^(1/(sig-1))/(1+(gt@{i})))*(exp(P@{i})/exp(P@{i}(+1)));
exp(chivar@{i}@{j}) = (rho_chi@{i}*exp(chiexo@{i}@{j}));
@#endfor
@#endfor


//Probability of adoption

@#for i in 1:n
@#for j in 1:n
exp(eps@{i}@{j})=epsbar@{i}@{j}*(exp(P@{i})*exp(Ha@{i}@{j})/exp(Yw))^betaa;
@#endfor
@#endfor


//Value of an innovation


@#for i in 1:n
exp(Vinnov@{i})=
@#for j in 1:n
   +exp(Jinnovator@{j}@{i})*exp(T@{i}(-1))/exp(T@{j}(-1))
 @#endfor
;
@#endfor


//FOC innovation    

@#for i in 1:n
exp(Hr@{i})=exp(P@{i})^(-1)*(exp(Vinnov@{i})*betar*lam@{i}*(exp(Yw))^(-betar))^(1/(1-betar));
@#endfor



//Law of motion of innovation

@#for i in 1:n
gz@{i}=lam@{i}*(exp(P@{i})*exp(Hr@{i})/exp(Yw))^betar*(exp(T@{i}(-1))/exp(Z@{i}(-1)));
@#endfor

//Law of motion of adoption

@#for i in 1:n
@#for j in 1:n
ga@{i}@{j}=exp(eps@{i}@{j})*(exp(Z@{j}(-1))*(1+gz@{j})/exp(A@{i}@{j}(-1))-1);
@#endfor
@#endfor


//Value of an unadopted technology

@#for i in 1:n
@#for j in 1:n
exp(J@{i}@{j})=-exp(P@{i})*exp(Ha@{i}@{j})*(exp(T@{i})/exp(A@{i}@{j}(-1)))*(exp(eps@{i}@{j})/(ga@{i}@{j}))+(1/exp(r))*(exp(eps@{i}@{j})*exp(V@{i}@{j}(+1))+(1-exp(eps@{i}@{j}))*exp(J@{i}@{j}(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt@{i})))*(exp(P@{i})/exp(P@{i}(+1)));
exp(Jinnovator@{i}@{j})=(1/exp(r))*(exp(eps@{i}@{j})*exp(Vinnovator@{i}@{j}(+1))+(1-exp(eps@{i}@{j}))*exp(Jinnovator@{i}@{j}(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt@{i})))*(exp(P@{i})/exp(P@{i}(+1)));
@#endfor
@#endfor




//FOC adoption

@#for i in 1:n
@#for j in 1:n
exp(P@{i})*exp(Ha@{i}@{j})*(exp(T@{i}(-1))/exp(A@{i}@{j}(-1)))*(exp(eps@{i}@{j})/(ga@{i}@{j}))=betaa*(1/exp(r))*exp(eps@{i}@{j})*(exp(V@{i}@{j}(+1))-exp(J@{i}@{j}(+1)))*((1+gt3)^(1/(sig-1))/(1+(gt@{i})))*(exp(P@{i})/exp(P@{i}(+1)));
@#endfor
@#endfor




//Total stock of knowledge

@#for i in 1:n
exp(T@{i})=
@#for j in 1:n
   +exp(A@{i}@{j})
 @#endfor
;
@#endfor


//Growth rates

@#for i in 1:n
(gp@{i})=(exp(P@{i}(+1))/exp(P@{i})-1)+((1/(1-sig))*gt3);
exp(TNT@{i})=gam@{i}*exp(T@{i});
@#endfor

@#for i in 1:n
(gz@{i})=(exp(Z@{i})/exp(Z@{i}(-1))-1)+gt3;
@#endfor

@#for i in 1:n-1
(gt@{i})=(exp(T@{i})/exp(T@{i}(-1))-1)+gt3;
@#endfor

(gt3)=exp(A31(-1))*(ga31)+exp(A32(-1))*(ga32)+exp(A33(-1))*(ga33);


@#for i in 1:n
(gtT@{i})=gt@{i};
(gtNT@{i})=gt@{i};
@#endfor



@#for i in 1:n
(gy@{i})=((exp(Y_L@{i}(+1)))/(exp(Y_L@{i}))-1)+(1+gt3)^(1/(sig-1))-1;
(gyy@{i})=((exp(Y@{i}(+1)))/(exp(Y@{i}))-1)+(1+gt3)^(1/(sig-1))-1;
@#endfor

@#for i in 1:n
@#for j in 1:n
(drp@{i}@{j})=(exp(rp@{i}@{j}(+1))/exp(rp@{i}@{j})-1)+(1+gt3)^(1/(sig-1))-1;
@#endfor
@#endfor



@#for i in 1:n
(gc@{i})=(exp(C@{i}(+1))/exp(C@{i})-1)+(1+gt3)^(1/(sig-1))-1;
@#endfor


@#for i in 1:n
@#for j in 1:n
(ga@{i}@{j})=(exp(A@{i}@{j})/exp(A@{i}@{j}(-1))-1)+gt3;
@#endfor
@#endfor

//Royalties
@#for i in 1:n
@#for j in 1:n
exp(rp@{i}@{j})=(exp(chivar@{i}@{j}))*exp(PiT@{i})*exp(A@{i}@{j}(-1))/exp(T@{i}(-1));
@#endfor
@#endfor

//Trade balance equation

@#for i in 1:n-1
exp(EX@{i})=
@#for j in 1:n
@#if j != i
+exp(pi_share@{j}@{i})*exp(P@{j})*exp(YT@{j})/(1+tau@{j}@{i}*(1-distexo@{j}@{i}))+exp(rp@{j}@{i})
@#endif
@#endfor
;
@#endfor


@#for i in 1:n-1
exp(IM@{i})=
@#for j in 1:n
@#if j != i
+exp(pi_share@{i}@{j})*exp(P@{i})*exp(YT@{i})/(1+tau@{i}@{j}*(1-distexo@{i}@{j}))+exp(rp@{i}@{j})
@#endif
@#endfor
;
@#endfor
@#for i in 1:n-1
exp(IM@{i})=exp(EX@{i})+B@{i}(-1)*exp(r(-1))*exp(Y@{i}(-1))*(1/exp(gy@{i}))-B@{i}*exp(Y@{i});
@#endfor

(w3)=(0);






//Productivity

 @#for i in 1:n
  exp(Y_L@{i})=exp(Y@{i});
@#endfor

 
end;


initval;

B1=0;
B2=0;
B3=0;

distexo11=0;
distexo12=0;
distexo13=0;

distexo21=0;
distexo22=0;
distexo23=0;

distexo31=0;
distexo32=0;
distexo33=0;


distexoNT11=0;
distexoNT12=0;
distexoNT13=0;

distexoNT21=0;
distexoNT22=0;
distexoNT23=0;

distexoNT31=0;
distexoNT32=0;
distexoNT33=0;



chiexo11=0.;
chiexo12=0.0;
chiexo13=0.0;

chiexo21=0.;
chiexo22=0.0;
chiexo23=0.0;


chiexo31=0.;
chiexo32=0.0;
chiexo33=0.0;



r=log((1+ss_g/(sig-1))/(1+ss_g/(1-sig))/beta);


C1=log(ss_c(1));
C2=log(ss_c(2));
C3=log(ss_c(3));

Y1=log(ss_y(1));
Y2=log(ss_y(2));
Y3=log(ss_y(3));

YT1=log(ss_yT(1));
YT2=log(ss_yT(2));
YT3=log(ss_yT(3));

YNT1=log(ss_yNT(1));
YNT2=log(ss_yNT(2));
YNT3=log(ss_yNT(3));




T1=log(ss_TT(1)/ss_TT(3));
T2=log(ss_TT(2)/ss_TT(3));
T3=log(ss_TT(3)/ss_TT(3));


TNT1=log(ss_TNT(1)/ss_TNT(3));
TNT2=log(ss_TNT(2)/ss_TNT(3));
TNT3=log(ss_TNT(3)/ss_TNT(3));


P1=log(ss_p(1));
P2=log(ss_p(2));
P3=log(ss_p(3));

PT1=log(ss_pT(1));
PT2=log(ss_pT(2));
PT3=log(ss_pT(3));


PNT1=log(ss_pNT(1));
PNT2=log(ss_pNT(2));
PNT3=log(ss_pNT(3));



EX1=log(ss_EX(1));
EX2=log(ss_EX(2));

IM1=log(ss_IM(1));
IM2=log(ss_IM(2));


Pi1=log(ss_pi(1));
Pi2=log(ss_pi(2));
Pi3=log(ss_pi(3));

PiT1=log(ss_piT(1));
PiT2=log(ss_piT(2));
PiT3=log(ss_piT(3));

PiNT1=log(ss_piNT(1));
PiNT2=log(ss_piNT(2));
PiNT3=log(ss_piNT(3));


w1=log(ss_w(1));
w2=log(ss_w(2));
w3=log(ss_w(3));

Vinnov1=log(ss_vinnov(1));
Vinnov2=log(ss_vinnov(2));
Vinnov3=log(ss_vinnov(3));

chivar11 = log(ss_chi(1,1));
chivar12 = log(ss_chi(1,2));
chivar13 = log(ss_chi(1,3));

chivar21 = log(ss_chi(2,1));
chivar22 = log(ss_chi(2,2));
chivar23 = log(ss_chi(2,3));

chivar31 = log(ss_chi(3,1));
chivar32 = log(ss_chi(3,2));
chivar33 = log(ss_chi(3,3));



Z1=log(ss_Z(1));
Z2=log(ss_Z(2));
Z3=log(ss_Z(3));

Y_L1=log(ss_YL(1));
Y_L2=log(ss_YL(2));
Y_L3=log(ss_YL(3));

Hr1=log(ss_hr(1));
Hr2=log(ss_hr(2));
Hr3=log(ss_hr(3));

gz1=ss_g;
gz2=ss_g;
gz3=ss_g;

gt1=ss_g;
gt2=ss_g;
gt3=ss_g;

gtT1=ss_g;
gtT2=ss_g;
gtT3=ss_g;

gtNT1=ss_g;
gtNT2=ss_g;
gtNT3=ss_g;


gy1=ss_g/(sig-1);
gy2=ss_g/(sig-1);
gy3=ss_g/(sig-1);

gc1=ss_g/(sig-1);
gc2=ss_g/(sig-1);
gc3=ss_g/(sig-1);

gp1=-ss_g/(sig-1);
gp2=-ss_g/(sig-1);
gp3=-ss_g/(sig-1);



ga11=ss_g;
ga12=ss_g;
ga13=ss_g;

ga21=ss_g;
ga22=ss_g;
ga23=ss_g;

ga31=ss_g;
ga32=ss_g;
ga33=ss_g;

Ha11=log(ss_ha(1,1));
Ha12=log(ss_ha(1,2));
Ha13=log(ss_ha(1,3));

Ha21=log(ss_ha(2,1));
Ha22=log(ss_ha(2,2));
Ha23=log(ss_ha(2,3));

Ha31=log(ss_ha(3,1));
Ha32=log(ss_ha(3,2));
Ha33=log(ss_ha(3,3));


rp11=log(ss_rp(1,1));
rp12=log(ss_rp(1,2));
rp13=log(ss_rp(1,3));

rp21=log(ss_rp(2,1));
rp22=log(ss_rp(2,2));
rp23=log(ss_rp(2,3));

rp31=log(ss_rp(3,1));
rp32=log(ss_rp(3,2));
rp33=log(ss_rp(3,3));


eps11=log(ss_eps(1,1));
eps12=log(ss_eps(1,2));
eps13=log(ss_eps(1,3));

eps21=log(ss_eps(2,1));
eps22=log(ss_eps(2,2));
eps23=log(ss_eps(2,3));

eps31=log(ss_eps(3,1));
eps32=log(ss_eps(3,2));
eps33=log(ss_eps(3,3));



V11=log(ss_v(1,1));
V12=log(ss_v(1,2));
V13=log(ss_v(1,3));

V21=log(ss_v(2,1));
V22=log(ss_v(2,2));
V23=log(ss_v(2,3));

V31=log(ss_v(3,1));
V32=log(ss_v(3,2));
V33=log(ss_v(3,3));


J11=log(ss_j(1,1));
J12=log(ss_j(1,2));
J13=log(ss_j(1,3));

J21=log(ss_j(2,1));
J22=log(ss_j(2,2));
J23=log(ss_j(2,3));

J31=log(ss_j(3,1));
J32=log(ss_j(3,2));
J33=log(ss_j(3,3));


Vinnovator11=log(ss_vinnovator(1,1));
Vinnovator12=log(ss_vinnovator(1,2));
Vinnovator13=log(ss_vinnovator(1,3));

Vinnovator21=log(ss_vinnovator(2,1));
Vinnovator22=log(ss_vinnovator(2,2));
Vinnovator23=log(ss_vinnovator(2,3));

Vinnovator31=log(ss_vinnovator(3,1));
Vinnovator32=log(ss_vinnovator(3,2));
Vinnovator33=log(ss_vinnovator(3,3));


Jinnovator11=log(ss_jinnovator(1,1));
Jinnovator12=log(ss_jinnovator(1,2));
Jinnovator13=log(ss_jinnovator(1,3));

Jinnovator21=log(ss_jinnovator(2,1));
Jinnovator22=log(ss_jinnovator(2,2));
Jinnovator23=log(ss_jinnovator(2,3));

Jinnovator31=log(ss_jinnovator(3,1));
Jinnovator32=log(ss_jinnovator(3,2));
Jinnovator33=log(ss_jinnovator(3,3));


pi_share11=log(ss_PishareT(1,1));
pi_share12=log(ss_PishareT(1,2));
pi_share13=log(ss_PishareT(1,3));

pi_share21=log(ss_PishareT(2,1));
pi_share22=log(ss_PishareT(2,2));
pi_share23=log(ss_PishareT(2,3));

pi_share31=log(ss_PishareT(3,1));
pi_share32=log(ss_PishareT(3,2));
pi_share33=log(ss_PishareT(3,3));



pi_shareNT11=log(ss_PishareNT(1,1));
pi_shareNT12=log(ss_PishareNT(1,2));
pi_shareNT13=log(ss_PishareNT(1,3));

pi_shareNT21=log(ss_PishareNT(2,1));
pi_shareNT22=log(ss_PishareNT(2,2));
pi_shareNT23=log(ss_PishareNT(2,3));

pi_shareNT31=log(ss_PishareNT(3,1));
pi_shareNT32=log(ss_PishareNT(3,2));
pi_shareNT33=log(ss_PishareNT(3,3));



A11=log(ss_AT(1,1)*exp(T1));
A12=log(ss_AT(1,2)*exp(T1));
A13=log(ss_AT(1,3)*exp(T1));

A21=log(ss_AT(2,1)*exp(T2));
A22=log(ss_AT(2,2)*exp(T2));
A23=log(ss_AT(2,3)*exp(T2));

A31=log(ss_AT(3,1)*exp(T3));
A32=log(ss_AT(3,2)*exp(T3));
A33=log(ss_AT(3,3)*exp(T3));

Yw=ss_Yw;
end;

steady;
check(qz_zero_threshold=1e-20);


load Final_values.mat;


endval;

rho_chi1 = ss_rhochi(1);
rho_chi2 = ss_rhochi(2);
rho_chi3 = ss_rhochi(3);

B1=0;
B2=0;
B3=0;

distexo11=0;
distexo12=0;
distexo13=0;

distexo21=0;
distexo22=0;
distexo23=0;

distexo31=0;
distexo32=0;
distexo33=0;


distexoNT11=0;
distexoNT12=0;
distexoNT13=0;

distexoNT21=0;
distexoNT22=0;
distexoNT23=0;

distexoNT31=0;
distexoNT32=0;
distexoNT33=0;



chiexo11=0.;
chiexo12=0.0;
chiexo13=0.0;

chiexo21=0.;
chiexo22=0.0;
chiexo23=0.0;


chiexo31=0.;
chiexo32=0.0;
chiexo33=0.0;



r=log((1+ss_g/(sig-1))/(1+ss_g/(1-sig))/beta);


C1=log(ss_c(1));
C2=log(ss_c(2));
C3=log(ss_c(3));

Y1=log(ss_y(1));
Y2=log(ss_y(2));
Y3=log(ss_y(3));

YT1=log(ss_yT(1));
YT2=log(ss_yT(2));
YT3=log(ss_yT(3));

YNT1=log(ss_yNT(1));
YNT2=log(ss_yNT(2));
YNT3=log(ss_yNT(3));




T1=log(ss_TT(1)/ss_TT(3));
T2=log(ss_TT(2)/ss_TT(3));
T3=log(ss_TT(3)/ss_TT(3));


TNT1=log(ss_TNT(1)/ss_TNT(3));
TNT2=log(ss_TNT(2)/ss_TNT(3));
TNT3=log(ss_TNT(3)/ss_TNT(3));


P1=log(ss_p(1));
P2=log(ss_p(2));
P3=log(ss_p(3));

PT1=log(ss_pT(1));
PT2=log(ss_pT(2));
PT3=log(ss_pT(3));


PNT1=log(ss_pNT(1));
PNT2=log(ss_pNT(2));
PNT3=log(ss_pNT(3));



EX1=log(ss_EX(1));
EX2=log(ss_EX(2));

IM1=log(ss_IM(1));
IM2=log(ss_IM(2));


Pi1=log(ss_pi(1));
Pi2=log(ss_pi(2));
Pi3=log(ss_pi(3));

PiT1=log(ss_piT(1));
PiT2=log(ss_piT(2));
PiT3=log(ss_piT(3));

PiNT1=log(ss_piNT(1));
PiNT2=log(ss_piNT(2));
PiNT3=log(ss_piNT(3));


w1=log(ss_w(1));
w2=log(ss_w(2));
w3=log(ss_w(3));

Vinnov1=log(ss_vinnov(1));
Vinnov2=log(ss_vinnov(2));
Vinnov3=log(ss_vinnov(3));

chivar11 = log(ss_chi(1,1));
chivar12 = log(ss_chi(1,2));
chivar13 = log(ss_chi(1,3));

chivar21 = log(ss_chi(2,1));
chivar22 = log(ss_chi(2,2));
chivar23 = log(ss_chi(2,3));

chivar31 = log(ss_chi(3,1));
chivar32 = log(ss_chi(3,2));
chivar33 = log(ss_chi(3,3));



Z1=log(ss_Z(1));
Z2=log(ss_Z(2));
Z3=log(ss_Z(3));

Y_L1=log(ss_YL(1));
Y_L2=log(ss_YL(2));
Y_L3=log(ss_YL(3));

Hr1=log(ss_hr(1));
Hr2=log(ss_hr(2));
Hr3=log(ss_hr(3));

gz1=ss_g;
gz2=ss_g;
gz3=ss_g;

gt1=ss_g;
gt2=ss_g;
gt3=ss_g;

gtT1=ss_g;
gtT2=ss_g;
gtT3=ss_g;

gtNT1=ss_g;
gtNT2=ss_g;
gtNT3=ss_g;


gy1=ss_g/(sig-1);
gy2=ss_g/(sig-1);
gy3=ss_g/(sig-1);

gc1=ss_g/(sig-1);
gc2=ss_g/(sig-1);
gc3=ss_g/(sig-1);

gp1=-ss_g/(sig-1);
gp2=-ss_g/(sig-1);
gp3=-ss_g/(sig-1);



ga11=ss_g;
ga12=ss_g;
ga13=ss_g;

ga21=ss_g;
ga22=ss_g;
ga23=ss_g;

ga31=ss_g;
ga32=ss_g;
ga33=ss_g;

Ha11=log(ss_ha(1,1));
Ha12=log(ss_ha(1,2));
Ha13=log(ss_ha(1,3));

Ha21=log(ss_ha(2,1));
Ha22=log(ss_ha(2,2));
Ha23=log(ss_ha(2,3));

Ha31=log(ss_ha(3,1));
Ha32=log(ss_ha(3,2));
Ha33=log(ss_ha(3,3));


rp11=log(ss_rp(1,1));
rp12=log(ss_rp(1,2));
rp13=log(ss_rp(1,3));

rp21=log(ss_rp(2,1));
rp22=log(ss_rp(2,2));
rp23=log(ss_rp(2,3));

rp31=log(ss_rp(3,1));
rp32=log(ss_rp(3,2));
rp33=log(ss_rp(3,3));


eps11=log(ss_eps(1,1));
eps12=log(ss_eps(1,2));
eps13=log(ss_eps(1,3));

eps21=log(ss_eps(2,1));
eps22=log(ss_eps(2,2));
eps23=log(ss_eps(2,3));

eps31=log(ss_eps(3,1));
eps32=log(ss_eps(3,2));
eps33=log(ss_eps(3,3));



V11=log(ss_v(1,1));
V12=log(ss_v(1,2));
V13=log(ss_v(1,3));

V21=log(ss_v(2,1));
V22=log(ss_v(2,2));
V23=log(ss_v(2,3));

V31=log(ss_v(3,1));
V32=log(ss_v(3,2));
V33=log(ss_v(3,3));


J11=log(ss_j(1,1));
J12=log(ss_j(1,2));
J13=log(ss_j(1,3));

J21=log(ss_j(2,1));
J22=log(ss_j(2,2));
J23=log(ss_j(2,3));

J31=log(ss_j(3,1));
J32=log(ss_j(3,2));
J33=log(ss_j(3,3));


Vinnovator11=log(ss_vinnovator(1,1));
Vinnovator12=log(ss_vinnovator(1,2));
Vinnovator13=log(ss_vinnovator(1,3));

Vinnovator21=log(ss_vinnovator(2,1));
Vinnovator22=log(ss_vinnovator(2,2));
Vinnovator23=log(ss_vinnovator(2,3));

Vinnovator31=log(ss_vinnovator(3,1));
Vinnovator32=log(ss_vinnovator(3,2));
Vinnovator33=log(ss_vinnovator(3,3));


Jinnovator11=log(ss_jinnovator(1,1));
Jinnovator12=log(ss_jinnovator(1,2));
Jinnovator13=log(ss_jinnovator(1,3));

Jinnovator21=log(ss_jinnovator(2,1));
Jinnovator22=log(ss_jinnovator(2,2));
Jinnovator23=log(ss_jinnovator(2,3));

Jinnovator31=log(ss_jinnovator(3,1));
Jinnovator32=log(ss_jinnovator(3,2));
Jinnovator33=log(ss_jinnovator(3,3));


pi_share11=log(ss_PishareT(1,1));
pi_share12=log(ss_PishareT(1,2));
pi_share13=log(ss_PishareT(1,3));

pi_share21=log(ss_PishareT(2,1));
pi_share22=log(ss_PishareT(2,2));
pi_share23=log(ss_PishareT(2,3));

pi_share31=log(ss_PishareT(3,1));
pi_share32=log(ss_PishareT(3,2));
pi_share33=log(ss_PishareT(3,3));



pi_shareNT11=log(ss_PishareNT(1,1));
pi_shareNT12=log(ss_PishareNT(1,2));
pi_shareNT13=log(ss_PishareNT(1,3));

pi_shareNT21=log(ss_PishareNT(2,1));
pi_shareNT22=log(ss_PishareNT(2,2));
pi_shareNT23=log(ss_PishareNT(2,3));

pi_shareNT31=log(ss_PishareNT(3,1));
pi_shareNT32=log(ss_PishareNT(3,2));
pi_shareNT33=log(ss_PishareNT(3,3));



A11=log(ss_AT(1,1)*exp(T1));
A12=log(ss_AT(1,2)*exp(T1));
A13=log(ss_AT(1,3)*exp(T1));

A21=log(ss_AT(2,1)*exp(T2));
A22=log(ss_AT(2,2)*exp(T2));
A23=log(ss_AT(2,3)*exp(T2));

A31=log(ss_AT(3,1)*exp(T3));
A32=log(ss_AT(3,2)*exp(T3));
A33=log(ss_AT(3,3)*exp(T3));

Yw=ss_Yw;
end;

steady;


perfect_foresight_setup(periods=1000);
perfect_foresight_solver;



shock13 = [0:0.1:1];
shockchi31= log([1:0.1:2.5]);//log([1:0.1:1.6]);
shockchi32= log([1:1:1]);//log([1:0.1:1.6]);
shockchi33= log([1:0.1:2.5]);//log([1:0.1:1.6]);
iterations = length(shock13)*length(shockchi31)*length(shockchi32)*length(shockchi33);

resultsC1 = zeros(iterations,1002);
resultsP1 = zeros(iterations,1002);
resultsgc1 = zeros(iterations,1002);
resultsC2 = zeros(iterations,1002);
resultsP2 = zeros(iterations,1002);
resultsgc2 = zeros(iterations,1002);
resultsC3 = zeros(iterations,1002);
resultsP3 = zeros(iterations,1002);
resultsgc3 = zeros(iterations,1002);

numit = 0;



for i=1:length(shock13)
    for j=1:length(shockchi31)
 for jj=1:length(shockchi32) 
for k=1:length(shockchi33)
                numit = numit + 1;                
                
                oo_.exo_simul(2:5,7) = [0];
                oo_.exo_simul(2:5,21) = [0];
                oo_.exo_simul(2:5,24) = [0];
                oo_.exo_simul(2:5,27) = [0]; 
                
                oo_.exo_simul(6:10,7) = shock13(i)/2;
                oo_.exo_simul(6:10,21) = shockchi31(j)/2;
                oo_.exo_simul(6:10,24) = shockchi32(jj)/2;
                oo_.exo_simul(6:10,27) = shockchi33(k)/2;

                oo_.exo_simul(11:1002,7) = shock13(i);
                oo_.exo_simul(11:1002,21) = shockchi31(j);
                oo_.exo_simul(11:1002,24) = shockchi32(jj);
                oo_.exo_simul(11:1002,27) = shockchi33(k);
             
                oo_ = perfect_foresight_solver_core(M_,options_,oo_);
                if oo_.deterministic_simulation.status
                    resultsC1(numit,:) = oo_.endo_simul(strmatch('C1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(7,:);
                    resultsP1(numit,:) = oo_.endo_simul(strmatch('P1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(14,:);
                    resultsgc1(numit,:) = oo_.endo_simul(strmatch('gc1',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(27,:);
                    resultsC2(numit,:) = oo_.endo_simul(strmatch('C2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(34,:);
                    resultsP2(numit,:) = oo_.endo_simul(strmatch('P2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(41,:);
                    resultsgc2(numit,:) = oo_.endo_simul(strmatch('gc2',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(54,:);
                    resultsC3(numit,:) = oo_.endo_simul(strmatch('C3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(61,:);
                    resultsP3(numit,:) = oo_.endo_simul(strmatch('P3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(68,:);
                    resultsgc3(numit,:) = oo_.endo_simul(strmatch('gc3',M_.endo_names,'exact'),1:1002);//oo_.endo_simul(81,:);
                    position(numit,:)=[numit i j jj k];
                end
end
end
            end
        end
